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Abstract 

The tempered evolution equation describes the trapped dynamics, widely appearing in nature, e.g., 
the motion of living particles in viscous liquid. This paper proposes the fast predictor-corrector ap¬ 
proach for the tempered fractional ordinary differential equations by digging out the potential ‘very’ 
short memory principle. The algorithms basing on the idea of equidistributing are detailedly described; 
their effectiveness and low computation cost, being linearly increasing with time t, are numerically 
demonstrated. 
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1 Introduction 

The tempered fractional calculus is a mathematical tool to describe the transition between normal and 
anomalous diffusions (or the anomalous motion in finite time or bounded physical space). The motion of 
a Levy flight particle can be characterized by the continuous time random walk (CTRW) model with a 
jump distribution function (f>{x) ^ < a < 2) and exponential waiting time distribution. The 

stable Levy measure for particle displacement makes arbitrarily large jumps possible and then results in 
divergent spatial moments M- However, the infinite spatial moments are not always appropriate for the 
physical processes; e.g., see Exponentially tempering the Levy measure is a popular way to 

make the moments of Levy distributions finite in transport models. Then we get the spatially tempered 
fractional Fokker-Planck equation mm- This paper focuses on the time tempered fractional derivative, 
which appears in the Fokker-Planck equation being derived from the CTRW model with tempered power 
law waiting time distribution US HSl HD]. The tempered power law waiting time measure has finite first 
moment and makes the trapped dynamics more physical; since sometimes it is necessary to make the first 
moment of the waiting time measure finite, e.g., the biological particles moving in viscous cytoplasm and 
displaying trapped dynamical behavior just have finite lifetime. The time tempered dynamics describes the 
phenomena of the coexistence/transition of subdiffusion and normal diffusion (or the subdiffusion in finite 
time) which was empirically confirmed in a number of systems |15L I16j . More applications for the tempered 
fractional derivatives and tempered differential equations can be found, for instance, in poroelasticity |12j . 
finance |2], ground water hydrology m, and geophysical flows [Ill- 

Tempered fractional calculus is the generalization of fractional calculus. Comparing with the classical 
differential equations, one of the big challenges we have to face when solving the (tempered) fractional 
equations is the expensiveness of its computation cost besides its complexity, this is because fractional 
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operators are pseudodifTerential operators which are non-local m- This paper focuses on providing the fast 
predictor-corrector approach for the following tempered fractional ordinary differential equation 


with the initial conditions 


0 


xit) = f{t,x{t)), 0<t<T, 


( 1 . 1 ) 


^ (e^*a;(t)) = Cfc, fc = 0,1, • • • , [a] - 1, 


( 1 . 2 ) 


where a G (0, oo), A > 0, Ck are arbitrary real numbers, and denotes the tempered fractional derivative 

in the Caputo sense [HIIS], defined by 




r, — Xt 


1 (i”(e^®u(s)) 


ds'^‘ 


ds, 


(1.3) 


T{n-a)Ja (t-s)“-"’+i 

where denotes the Caputo fractional derivative m Note that when A = 0, the Caputo tempered 
fractional derivative reduces to the Caputo fractional derivative. 


One of the effective and popular methods for numerically solving the initial value problem (1.1)-(1.2) with 


A = 0 is the predictor-corrector approach, raised by Diethelm et al mw- This method is improved in the 
Preliminary and Appendix sections of [S] (or see the review article 0), where around half of the computation 
cost seems to be cut off and the convergence order is improved from min{l -|- a, 2} to min{l -|- 2 q;, 2}. Some 
efforts have been made to reduce the computation cost by using the nested meshes basing on the fixed memory 


principle [H] and the short memory principle [10] of fractional operator when a G (0,1) in ( |1.1| ), and the 
short memory principle is apprehended from a new point of view in [3] which extends its effective range to 
a G (0, 2). With the nested meshes, the computation cost can be reduced from 0{h~'^) to 0{h~^ log{h~^)) 
while not losing the numerical accuracy, but its seems that this is more theoretical claims rather than 
numerical practices; in fact the numerical simulation results in [4l[T0] show this. More recently, the so called 
Jacobian-predictor-corrector approach, is introduced by Zhao and Deng [22], in which the accuracy is greatly 
improved while the computation cost is not increased. 

In this paper, we dig out the potential of the short memory principle of tempered fractional operators 
and apply it felicitously to reduce the computation cost by using the idea of equidistributing meshes [21] for 
numerically solving the initial value problem ([LT|)-(|r^ with the predictor-corrector method. In the next 
section, we introduce the techniques of equidistributing meshes, present the detailed numerical schemes, and 
describe the algorithms by pseudo codes. Numerical examples to illustrate the efficacy of the algorithm are 
given in Section]^ We conclude the paper with some remarks in the last section. 


2 Predictor-corrector Algorithms with Equidistributing Meshes 


The initial value problem (1.1)-(^1.2| is equivalent to the following Volterra integral equation 


M-i -xt 




g-Af^fc 1 




fc =0 


fc! r(a) 


{t-rT /(T,a;(r))dr 


:= + [ e V('r, a;(r))dr, 

r(a) Jo 


( 2 . 1 ) 


M-i 


where xo{t) := ■ For uniform nodes tn+i = {n + l)h, n = 0,1, - ■ ■ ,N, with h = T/N being 

k—0 

the steplength of numerical computation, the above equation can be recast as @10 

x(tr,+i) = xo(tn-Hi) + - r)“"V(Da;(T))(iT 

r(a) Jt„ 




( 2 . 2 ) 


-A(t„ + 1-T) 


(<„+i-r)“ f{T,x{T))dT, 
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or 


X{tn+l) = X{tn) + Xo{t„+l) - Xo{t„) 

1 r' 


+ 


e ^\tn+l-T)°‘ ^f{T,x(T))dT 


T{a)J,^ - (2.3) 

+ fR C - r)“-i) f{T,x{T))dT. 

For the single step integral (t„+i — T)““^g(T)dr, we can employ the rectangle or trapezoidal 

quadrature formula to approximate it, i.e., 


e-^(‘’'+i-")(t„+i - r)“-i5(T)dT 


or 


/5+^(t„+i - r)“-i g(«^+i)("-*^)+«~"S(«^)d^+i-") dr 
[ae-^''g(t„) +5(t„+i)] . 


(2.4) 


(2.5) 


So, obviously, for the numerical approximations of (2.2) and (2.3), the main computational cost comes from 


the term - dr. Luckily when t„+i —)■ +oo, (t„+i — r)““^ decays algebraically with the order 1 — a 
for a S (0,1], and (t„+i — t)““^ — (t„ — t)““^ decays with the order 2 — a for a G (0,2), these are 
the so-called short memory principle m and the one apprehended from a new point of view both 
the integral kernels — t)““^) and — t)““^ — — r)““^) ex¬ 

ponentially decay. We felicitously apply these decay properties to reduce the computation cost based on 
equidistributing meshes rather than roughly using the nested meshes. The linearly increasing computation 
cost with time t is shown in the following simulations. 


2.1 Algorithms for (2.2) when a G (0,1] 


In this subsection, we design numerical schemes for a G (0,1] by using (2.2) with the kernel e {tn+i — 


As mentioned above, how to compute the integral — x{T))dT efficiently 

is the key of reducing the computation cost to obtain the numerical approximations to x{tn+i)- In the 
previous predictor-corrector work introduced in laiHiE], as well as its improved version in [HE], n steps’ 
calculations are used to approximate the integral f*" - dr. While, by noticing that (t„+i — t)““^ decays 
with power 1 — a, and g-Htn+i-r) exponentially (see figures below), we can actually select less mesh 

points, say 


0 — '^0,n ^ '^l,n ^ * * * *\ Xjj 


— It, 


( 2 . 6 ) 


at which to approximate the term • dr by compounding two-point trapezoidal quadrature formula. In 
the following, we list two ways of choosing the mesh points 


2.1.1 Two equidistributing options for selecting the quadrature nodes 

Assuming that we have already get the points Ti.„, 0 < i < to„. For not losing the accuracy but reducing 
the computation cost, the first way of selecting the next point is based on the principle that the values 

of the function ?/(t) = {tn+i — r)““^ are equally distributed, i.e.. 


y{Ti+l^n) - y{Ti,n) 

= g-A(t„+i-n,„-(ri+i,„-n,„)) (t„+i 

= Ay, 




O' — 


1 


(2.7) 
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where Ay is a given small positive real number. Since (2.7) is a nonlinear equation w.r.t. we use its 

linear part instead to select the wanted point, and add another condition that promise is at least one 

step away from r^. That is, let 


= max • 


solve(fj+i,„ - n^n = h, n+i^n), 1 

solve((fj+i,„ - n^n) [(1 - a) + A(t„+i - n^n)] = J ’ 


( 2 . 8 ) 

where ‘solve(equ, var)’ means the solution of ‘equ’ with unknown variable ‘var’. Secondly, to avoid involving 
non-equal divided nodes, we take 






• h. 


(2.9) 


Figs. A and Algorithm illustrate this equal-height distribution method more clearly. 



Figure 2.1. The relationship between the selected nodes {rt} and the function y = e -’'(*"+1 '^)(t„_|_i — t)“ ^ 
in the algorithm of equal-height distribution, where Ay = h = 1/10, = 2, A = 1, and a = 0.2. 

The second way of choosing the mesh points {n^n} is to make the integrations of y(r) = g-Htn+i-r) 
w.r.t. T on any interval are almost the same, i.e., 


m + l.r, 


„-A(t„ + i-T) 


(t„+i - t)“ ^dr = As, 


( 2 . 10 ) 


where As is a given small positive real number. Also, since (2.10) is a nonlinear function w.r.t. we 

simply approximate it by 


rn+i,, 


2-^(‘"+^-")(tn+i-r)“-idr 


{tn+l - 






( 2 . 11 ) 


Thus, 


n+i,n = max • 


SOlve(fi+p„ - Ti^n = h,fi+i^n), 
solve ((fj+i,„ - Ti,„) = As(t„+1 - 


( 2 . 12 ) 
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Figure 2.2. The relationship between the selected nodes {rt} and the function y = e — t)“ ^ 

in the algorithm of equal-height distribution, where Ay = h = 1/10, = 2, A = 1, and a = 0.8. 


and then by taking 


”^2+1,n 


"^ 2 + 1,71 

h 


■ h, 


(2.13) 


Ti^i^n belongs to the set of the uniform grid points 
equal-area distribution criterion of selecting the mesh point 


Figs. 2.3||2.4 and Algorithm 


show this 


} more concretely. 


Remark 1. From Figs. 2. 1^2.4 we can see that since the function yi^r) = —steepens 

with the decrease of a, for fixed h and fixed Ay in the equal-height distribution method or fixed As in the 
equal-area one, the less a is, the fewer the quadrature nodes {xi^n} axe selected. This is different from the 
predictor-corrector approach proposed in 00[I or the improved version in m, in which the integral /q" -dr 
is always approximated on all the mesh points nodes number is unrelated to the value of a. 


2.1.2 Predictor-corrector approach for ( |2.2[ ) when a G (0,1] 

Since the mesh nodes {xi^n} chosen from Algorithm or Algorithm still belong to the set of the uniform 
nodes, we denote 

Xi,n=tni, z = 0,1, • • • ,m„, (2.14) 
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Figure 2.3. The relationship between the selected nodes {rt} and the function y = e -'‘(*"+1 — t)“ ^ 

in the algorithm of equal-area distribution, where h = 1/10, As = h, tn+i = 2, A = 1, and a = 0.2. 



Figure 2.4. The relationship between the selected nodes {rt} and the function y = e -’'(*"+1 — t)“ ^ 

in the algorithm of equal-area distribution, where h = 1/10, As = h, tn+i = 2, A = 1, and a = 0.8. 
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Algorithm 1 The pseudocode of equal-height distribution algorithm for (2.8)-(2.9) 


* = 0 ; Ti^n = 0 ; % ro,n = 0 ; 

Tc = 0 ; % current node is To,n 


while Tc < t 


if 


do 

[(l-a)+A(t„+i-rc)] ' 

> tn then 


break; 

end if 


% if y = — t)““^ changes too fast 

if Ti+i^n — Tc < h then 

Ti+i,n = Tc + h; % make Ti+i^n be one step away from Ti^n 

else 


Ti+i,n = \ji+i,nlh\ * h] % let belong to the uniform mesh points {tj}f^Q 

end if 

'^C - I'z-t-1,715 i - "2 “t“ 1, 

end while 


where tj = jh. Note that = 0 and t^^^ = t^. In this way, by using the product trapezoidal quadrature 
formula to replace the integral, there is 

f (<„+i - T)“"^g(T)dr 


r'T'i+lp 




2=0 ' 
rrin — l 

E 

i=o 

7n„-l pt 


,a_i e-^(*"+i-^-+i-)5('r7+i,r7)(r - t,,„) -f - r) , 

(^n+l '^) 

i n '^i,n 


= E 


(t„+i - t) 


»-ie -f e ^(‘"+1 g{tn,)itni+, - t) 


2=0 

m ^-1 X{tr^+l-tr^ 


= E 




2=0 


+ E 


2=0 


(ni+i-ni)h j 

g-A(t„+i-t„.)g(^^J rt 

(rii+i - ni)h 


+ l ty, 




(t„+i-r)“ ^(t„ -T)dr 


dr 


= E 


/i“e-A(*-+i-‘".)g(t„J 


2=1 

rrin — 1 


+ E 

2=0 

rrin — l 

= E 


(Uj - n7-i)(n -I- 1 - _ (n -I- 1 - rty)“+^ - (n -|- 1 - ni-i)°+^ 

a Q;(a -f 1 ) 

{uj+i - ni)(n -I-1 - ni)°‘ ^ (n -|- 1 - ny +i)“+^ - (n -|- 1 - ni)“+^ 


{Ui -Ui-i) 

(rii+i - Ui) [ a ' a(a -f 1) 

a(a + 1) 

(n -I- 1 — — (n -I- 1 — (n -I- 1 — — (n -f 1 — ni_i)“+^ 


^ 2+1 ^2 


^2 ^ 2—1 


+ 


h“e 

nia(a -I- 1) 


{(n -I- 1 - ni)“+^ - (n -I- l)“[n -f 1 - (a -I- l)ni]} 


+ 7 - {(n -f 1 - nyyy„_i)“+^ ” (o "f !)(« - nm„-i) - 1 } . 

[n — njn^_i)a[a + 1) 


(2.15) 
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Algorithm 2 The pseudocode of equal-area distribution algorithm for (2.12)-(2.13) 


* = 0 ; Ti^n = 0 ; % ro,n = 0 ; 

Tc = 0; % current node is To,n 

while Tc < tn do 

n+i^n =Tc + As(t „+1 - rc)'^-agHt„+i-T,). 
if n+i,n > tn then 


break; 

end if 


% if y = — t)““^ changes too fast 

if Ti+i^n — Tc < h then 

Ti+i,n = Tc + h; % make Ti+i^n be one step away from Ti^n 

else 


Ti+i,n = \ji+i,nlh\ * h] % let belong to the uniform mesh points {tj}f^Q 

end if 

'^c — "^z-t-l.ni i — "2 “t“ 1, 

end while 


Assume that the approximations Xj « x(tj), j = 1,2, ■■■ ,n, have been obtained. Then by combining with 
the one step approximations (2.4) and (2.5), the predictor-corrector approach to compute Xn+i ~ x{tn+i) 
for a G (0,1] can be yielded as 
Case 1 (n = 0): 


ccf’’ = a;o(ti) -f 
xi = xoih) + 


e-’'^f{0,xo), 


r{a + l) 

h°‘ 


(2.16) 


Case 2 (n > 1): 


= Xo{tn+l) + 


^n+1 — ^o(^n+l) T 


r(a + 2 ) 

h°‘ 

r(a-H 2 ) 

/i“ 

r(a-h 2 ) 


[f{h,x^'^)+ae ^'*/( 0 ,xo)]; 


^ ^ ^i,n-\-lfiti,Xni')~\~bnf(tn,Xn^ 

. i=0 

^ 0‘i,n+lf {tj, Xm) T f {tn+li Xn^i) 


(2.17) 


.i=0 


where 


-X(„ + l)h 


{(n -I- 1 — ni)“+^ — {n + l)“[n -|- 1 — (a -I- l)ui]} , 


^z,n+l — ^ 


and 




(n+1 —— (n+1 —Tii 

rii+i —rii 


(n+1 ——(n+1—ni i)° 


z = 0, 


1 < i < rrin — 1 , 


- [(n -I- 1 - nm„-i)“+^ - (n -I- 1 - nm„-i)] 


i = TO„; 


— ^rrin ,n-\-l 




(2.18) 


(2.19) 


If no equidistributing principle is adopted, i.e., Ui = i, the predictor-corrector approach for ( |2.2[ ) to 
compute Xn+i ~ x{tn+i) is described as 

Case 1 (n = 0): being the same as (2.16); 































Case 2 (n > 1): 


=xo{tn+i) + 

Xn-\-l — ^o(^n+l) “t“ 


r(a + 2) 

r(a + 2) 


d,^n+lf{ti, Xi) + e - 1) • f(tn, Xn) 

.i=0 

n 

^ dj^n+lf {tjj Xj) f [tn+l, X^_^_l) 


Li=0 


where 


di^n+l — 


e ''‘("+1)^ — (n + l)“(n — a)) , i = 0, 

e-^("+i-*)'* • ((n - i)“+i - 2(n + 1 - i)“+i + (n + 2 - i)“+i) , 1 < i < n. 


( 2 . 20 ) 


( 2 . 21 ) 


Theorem 1. For the tempered fractional initial value problem assume f(t,xit)) of belongs 

to for some suitable T. Then for the scheme l2.20)-(2.21 ) we have 


max \x{tn)-Xn\ = 
0<n<N 


OQif), if C( > 0.5, 
0(/ii+2“), if 0<a< 0.5. 


Proof. Combining Appendix A of [5] with Theorem 2.5 and Lemma 3.1 of [9], we can obtain the result of 
this theorem. □ 


2.2 Algorithms for (2.3) when a G (0,2) 


Now we apply the predictor-corrector algorithms to Eq. (2.3) for a € (0,2). To begin, we analyze the 
properties of the kernel function y(T) = (<„+i — r)““^ — e~^^*"'~'^\tn — t)““^. On one hand, by 

noticing A > 0, there is 


< e-^O.+1-O . (2.22) 

together with the short memory principle [13 E] that (t„+i — t)““^ — (t„ — r)““^ decays with the order 
2 — a, we can see that ?/(r) decays faster than exponential order A for fixed r when —>■ cx). On the other 
hand, it can be seen from 

' = e-^^x^-^ [(a - 1) - Ax] (2.23) 

that ?/(t) < 0 if a S (0,1 ]; and if a € ( 1 , 2 ) and when tn+i — t< i.e., t > there is ?/(t) > 0; 

while if a G (1, 2 ) and when t„ — r > i.e., r < t„ — there is ?/(r) < 0. Also, 2 /(t) is lower bounded. 

Specifically, 


y{r) = e-^(‘"+^-")(t„+i-r)“-i-e-^(‘"-^ 


> {tn - tY 


g-A(t„+i-r) _ g-A(i„-T) 


(tn - tY 


= (t„-r)“-ie-^(‘"+i-")(l-e^'‘) 

= -(tn-r)“-ie-^(‘"+^-") • 0 (h), 


(2.24) 


that is, 2 /(t) is very close to zeros when it is negative; thus, we can directly neglect the part when j/(t) < 0 . 

Basing on this fact, similar to the case of dealing with \2.2\ when a G (0,1], we can explore more effi¬ 
cient ways to approximate the integral Jp" — t)““^ — '^\tn — t)°‘ g(T)dT with 

trapezoidal quadrature formula. Corresponding to Algorithm [^ and Algorithm we provide equal-height 
as well as equal-area distribution methods for choosing the quadrature nodes in the subsection below. 
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2.2.1 Two equidistributing options for selecting the quadrature nodes 

Suppose that those to be chosen points are sequential as in ( |2.6[ ), and we have already got the points 

Ti.n, 0 < i < TO„. By the analysis above, we further let > min{t„,t„+i — for a € (1,2). Like the 
equal-height distribution principle in Subsection 2.1, an intuitive way of selecting the next point is, 

firstly, let 


„-A(t„+i—fi 




(tn+1 


+ (-1)^“^ Ay, 


(2.25) 


where Ay is a given small positive real number, and if a G (0,1), ~Ay is used, else, if a G (1,2), it is 
replaced by Ay. That is. 


A[t„+1 Ti,„ (Ti+I,,, Ti,„)] 

g-A[t„-n,„-(n+i,„-ri,„)] _ 


-n,n)]“ ^-e 1 

Tj.„)]““^ - _ ri,„)““i -f (-l)rfo Ay. 


(2.26) 

Since the above equation is nonlinear w.r.t. we use its linear part instead to select the wanted point. 

That is. 


Ti+i^n = max 


solve '^i,n — h, T^-|-l^7^) , 

solve((fi+q„ - Ti^n) [A(t„+i - 'ri_„)““i - (a - l)(t„+i - rj,„)““^] 
= (fi+i,„ - Ti^n)e^^ [A(t„ - Ti,„)““^ - (a - l)(t„ - Ti,„)““2] 

+ (_l)folAye^(‘.+i-n,„)^^^^^_^) 


(2.27) 


Then, to avoid being non-equal divided nodes, we take 

^z+l,r 




• h. 


(2.28) 


The above descriptions are demonstrated more clearly by Figs. |2.5||2.8| and Algorithmic 

Similar to the equal-area idea in Subsection 2.1, here for the second way of choosing the mesh points 
{Ti.ra}, we expect firstly that 


[e-^(‘"+^-")(t„+i - r)“-i - e-^(*"-")(t„ - r)“-i] dr = (-1)^“^ As, 


(2.29) 


where As is a given small positive real number, and if a G (0,1), —As is used, else, if a G (1, 2), it is replaced 
by As. Again, since it is a nonlinear equations w.r.t. we approximate it by 


(^n-t-l 7t,n) 


+ 

_ 


/ Ti 

■ 


dr = (—1)1^“^ As. 


So, let 


solve (fj+i,„ - n^n = h, n+i^ri) , 

fi+i^n = max <1 solve(A(fi+i,„ - n,„) [(t„+i - ri,„)““i - - Ti,„)““C 

= (-1) foT AsAe^(*"+i"'^-."), 


Then by taking 




'^2 + 1,7 


(2.30) 


(2.31) 


(2.32) 


n belongs to uniform grid points have a better view of this equal-area distribution 

criterion of selecting the mesh from Figs. 2.9|2.12 and Algorithm |C 
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Figure 2.5. The relationship between the selected nodes {ri} and the function y = (t„ — — 

in the algorithm of equal-height distribution, where h = 1/10, Ay = h/5, t„+i = 2, 

A = 1, and a = 0.2. 



Figure 2.6. The relationship between the selected nodes {ri} and the function y = — t)““^ — 

_r)““i in the algorithm of equal-height distribution, where h — 1/10, Ay = /i/5, t„+i = 2, 
A = 1, and a = 0.8. 
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Figure 2.7. The relationship between the selected nodes {Ti} and the function y = (t„+i — 

— r)““^ in the algorithm of equal-height distribution, where h = 1/10, Ay = h/5, 
tn+i = 5, A = 1, and a = 1.2. 



Figure 2.8. The relationship between the selected nodes {ri} and the function y = (t„+i — 

— t)““^ in the algorithm of equal-height distribution, where h = 1/10, Ay = h/5, 
tn+i = 5, A = 1, and a = 1.8. 
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Figure 2.9. The relationship between the selected nodes {ri} and the function y = (t„ — — 

— r)““^ in the algorithm of equal-area distribution, where h = 1/10, As = h/5, tn+i = 2, 
A = 1, and a = 0.2. 



Figure 2.10. The relationship between the selected nodes {ri} and the function y = — t)““^ — 

algorithm of equal-area distribution, where h = 1/10, As = h/5, tn+i = 2, 

A = 1, and a = 0.8. 
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Figure 2.11. The relationship between the selected nodes {r^} and the function y = (tn+i — 

— r)““^ in the algorithm of equal-area distribution, where h = 1/10, As = h/5, 
tn+i = 5, A = 1, and a = 1.2. 



Figure 2.12. The relationship between the selected nodes {t^} and the function y = (t„+i — 

— t)““^ in the algorithm of equal-area distribution, where h = 1/10, As = h/5, 
tn+i = 5, A = 1, and a = 1.8. 
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Algorithm 3 The pseudocode of equal-height distribution algorithm for (2.27)-(2.28) 


* = 0; Ti^n = 0; % ro,n = 0; 

Tc = 0; % current node is To,n 

if 0 < a < 1 then 

Ay = — Ay; % if a G (0,1], take negative value 
else 


Ti.n = tn +1 — % start from the point where y{T) > 0 

Ti,n = ^ belong to the uniform mesh 

i = 1 ; Tc = Ti_„; % current node is 

end if 

while Tc < tn do 

n+i^n -TcG- [A(t„+i-r,) —i-(c<-l)(t„+i-r,) —2]-e^''[A(t„-r,)“-i-(a-l)(t„-r,)“-2] > 

if T^+i,n > tn then 

break; 

end if 

% if y = — t)““^ — e~^^^^~'^\tn — r)““^ changes too fast 

if Ti+i „ — Tc <h then 

Ti+i,n = Tc + h; % make Ti+i^n be one step away from Ti^n 

else 


Ti+i,n = [Ti+i^n/h\ * h] % let Ti+i^n belong to the uniform mesh points {tj}”=o 

end if 

'Cc — i — "2 “t“ 1 , 

end while 


2.2.2 Predictor-corrector approach for ( |2.3[ ) when a G (0, 2) 

We still denote the mesh nodes chosen from Algorithmor Algorithmas done in (2.14). Also, by 

noticing that tn^ = 0 and from (2.16), we can get 


e-A(‘"+-")(t„+i - t)“-i - g{T)di 


= E 

rrin — 

- E 


m+lo 




e ^g{T)dT- ^ / e ^\tn - t^ ^g{T)dT 


i=0 •'T 
rrin —1 

- E 

2 = 0 ' 


rt,n i^Q ^Ti^ri 

(tyj+l "T” j (^'7' 

ri n ^2+1,n "^2,72 


^2+1,n ^2,r; 


7 I 7T / . (^i.n-l-l ~ Ci n+l)<7(tni)) 

(a + 1 ) 

^ ’ i=0 


(2.33) 


where {ai^„+i}(EQ are defined in (2.18), and the definitions of Ci^n+i, i = 0 , 1 ,--- ,to„, are obtained by 
replacing n + 1 of (2.18) corresponding to f = 0,1, • • • , m„, with n. However, it should be noted that 

^ 2 , 71+1 7 ^ because 7 ^ (rz l)fc for k — 1 , 2 , * * *. 

Assume that we have calculated the approximations Xj « x{tj), _) = 1, 2, • • • , rz. Then by combining with 
the one step approximations (2.4) and (2.5), the predictor-corrector approach to compute Xn+i ~ a:(tn+i) 
for a G (0, 2) can be yielded as 
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Algorithm 4 The pseudocode of equal-area distribution algorithm for (2.31)-(2.32) 


* = 0; = 0; % ro,n = 0; 

Tc = 0 ; % current node is To,n 

if 0 < a < 1 then 

As = — As; % if a G (0,1], take negative value 
else 


Ti^n = tn +1 — % start from the point where yfr) > 0 

Ti,n = \ji,nlh\ ■ h] % let belong to the uniform mesh 
i = 1 ; Tc = Ti_„; % current node is 

end if 

while Tc < tn do 

_I 

Ti+l,n —Tc-r [(t„+i-Te)“-i-e^'*(t„-T,)“-i] ’ 

if n+i,n > tn then 

break; 

end if 

% if y = — t)““^ — e~^^^^~'^\tn — r)““^ changes too fast 

if Ti+i „ — Tc <h then 

Ti+i,n = Tc + h; % make Ti+i^n be one step away from Ti^n 

else 


Ti+i,n = \ji+i,nlh\ * h] % let Ti+i^n belong to the uniform mesh points {tj}”=o 

end if 

'^c — "^z-t-l.ni i — "2 “t“ 1 , 

end while 


Case 1 (n = 0): 


= Xoih) + p. \ ,x e ^^f{0,xo), 

L[a + l) 

< (2.34) 

xi = xo{ti) + ^ [f{h, xf’') -f ae"^^/( 0 , xq)] ; 

i (a -f 2) 


Case 2 (n > 1): 


y,Pr 


= Xo(tn+l) -\-Xn- Xo{tn) + 


r(a-h2) 

rrin — l 

f (j'i i Xm ) “1“ ,n+l If (fn i Xn) 

, 2=0 


Xn-\-l — ^o(^n+l) H“ Xn XQ(tn) H“ 


r(a + 2 ) 

m„ — 1 

) “ 1 “ , 21+1 ^22171 ,n+l“l“^^ (^n; ^n) “ 1 “ f (^n+1 5 ^n+ 1 ) 

L 2=0 


(2.35) 


3 Numerical Examples 


Example 1. Our first example deals with the case that the unknown solution x(t) has a smooth derivative 
of order a. Specifically, we shall consider the following equation as in HU- 


C 

0 


^a.A 


x{t) 


= e 


-xt 


r(9) s-c _ r(5 + a/2) ^4_„/, 
r(9-a) r(5-a/2) 


^r(a + l) + (V2-t4)^-(e^+)3/^ 


(3.1) 
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with the initial condition(s) x( 0 ) = 0 (and \e^'*‘x(t)\' = 0 if 1 < a < 2). 

The exact solution of this initial value problem is 

x{t) = e"^* , 


(3.2) 


so, the right function ^Df’^x^t) = t^~°‘ — + |r(Q; + 1) € C'^[0, T] for arbitrary T > 0 

and a G ( 0 , 2 ). 

Tablej^verifies that, by choosing suitable parameters Ay or As for different a G (0,1], the equidistributing 
predictor-corrector methods can reach to the convergence order minjl -|- 2a, 2} as given in Theorem 1. Also, 
it can be seen that when the time interval is not very large (T = 1) and A = 1, the equal-area distribution 
predictor-corrector methods have already shown their benefits in computation cost for a G (0,1]. As for the 
case of a G (1,2), both equidistributing methods show their huge superiority in the aspects of numerical 
accuracy as well as computation cost. Noting that to compute Xn+i, the least three required quadrature 
nodes are 0, tn, and tn+i, except when n = 0, only two points 0 and ti are employed. So, given h and T, the 
total times of the quadrature nodes being used in the procedure is 2 -|- 3(^ — 1), i.e., 29, 59, 119, 239, and 
479, respectively, when h = 1/10, 1/20, 1/40, 1/80, 1/160, and T = 1, which coincides with the values of M 
in Table while the total times of the quadrature nodes being used for nonequidistributing scheme (2.161 
and ( 2.20\ - ^2(2i ] are 65, 230, 860, 3320, and 13040, separately. That means, when a G (1,2), while keeping 
high numerical accuracy, the computation cost of the equidistributing methods are linearly increasing with 
time, comparing with the 0{h~^) expenditure in the predictor-corrector methods of [7l[5]. 

Remark 2. For not losing accuracy, the basic strategy of choosing the values of Ay or As is: 1) for the 


algorithms generated from (2.2) where a G (0,1), if a is close to 1, Ay or As is approximately equal to h, 
and if a is close to 0, Ay or As can be bigger than h, say, lO/i; 2) for the algorithms generated from l[2.t^ 
where a G (0,2), when a is closer to 1, Ay or As should be smaller (much smaller than h if a G (0,1) and 
around h if a G (1,2)/. 


Example 2. Next we come to the case that the unknown solution x{t) itself is a smooth function, but the 
given function f(t,x(t)) has weak regularity. Specifically, we consider the linear equation: 


C 


DTx{t) = 


„-\t 




r(3-c 

2 


^2-a _ -t 


^2-a _ 


t^ a _ gAtj-j-^^ _|_ ^2 _ ^ 


r(3-a)‘' r(2-a)' 

with the initial condition(s) x( 0 ) = 0 (and [e^^x{t)\' = —1 z /1 < a < 2 /. 

The exact solution of this initial value problem is 


for a > 1, 
for a < 1, 


(3.3) 


x{t) = e ^*{t^ — t), 


(3.4) 


so, the right function 


U i~)a 
0 


x(t) = 


2 g2-a 

r(3-Q)'' ’ 

2 +2-a 

r(3-Q)'' 


1 j-l-g 

r( 2 -a) '' 


for a > 1 , 
for a < 1 , 


(3.5) 


is continuous, but its first order derivative is infinite at t = 0 . 

Table shows that when the time interval is larger (T = 5), the advantage of computation cost for 
equidistributing schemes becomes more obvious. For a G (1, 2), the proposed methods of this paper possess 
more obvious advantages; they can even reach to the convergence order of 2 , although the function f(t,x(t)) 
has a weak regularity. Moreover, as discussed in Example the total times of the quadrature nodes being 
used are nearly equal to 2 + ‘i{'^ — 1 ) for a G ( 1 , 2 ), which means the computation cost of the equidistributing 
methods increases linearly with the time evolution. We can have a better view from Fig. |3.1| with a = 0.5 
and Fig. |3.2| with a = 1.5, where the equidistributing methods manage to linear growth of CPU time for 
long time interval (T = 50). In this section, all numerical computations are done in Matlab 7.11.0 on a 
normal laptop with 1GB of memory. 
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Table 1. The maximum errors {cmax), convergence rates (CO), and the total times of quadrature nodes 
being used (M) of Example at T = 1 and A = 1, by using the equal-height distribution as well as the 
equal-area distribution ones, for different a, h, Ay, and As, respectively. 


a 


a — 0.2 (Ay 

' = As 

= lOh) 

a = 0.5 (Ay = As 

_ 5/iA 

2 ' 

a = 0.8 (Ay = 

As = h) 

Method 

h 

^max 

CO 

M 

^max 

CO 

M 

^max 

CO 

M 

Equal- 

1/10 

5.83 le-1 

- 

32 

2.87 le-2 

- 

57 

8.23 le-3 

- 

64 

height 

1/20 

2.01 le-1 

1.52 

148 

6.29 le-3 

2.19 

205 

1.62 le-3 

2.35 

228 


1/40 

4.04 le-2 

2.32 

588 

1.44 le-3 

2.12 

777 

3.44 le-4 

2.23 

854 

of (2.2) 

1/80 

9.59 le-3 

2.08 

2327 

3.42 le-4 

2.08 

3034 

8.22 le-5 

2.07 

3300 


1/160 

2.78 le-3 

1.78 

9214 

8.24 le-5 

2.05 

11938 

1.98 le-5 

2.06 

12976 

Equal- 

1/10 

5.84 le-1 

- 

29 

3.29 le-2 

- 

43 

8.23 le-3 

- 

61 

area 

1/20 

2.00 le-1 

1.55 

59 

6.80 le-3 

2.28 

160 

1.62 le-3 

2.35 

221 


1/40 

4.07 le-2 

2.30 

287 

1.56 le-3 

2.12 

609 

3.34 le-4 

2.27 

830 

of (2.2) 

1/80 

9.68 le-3 

2.07 

1206 

3.64 le-4 

2.10 

2379 

7.23 le-5 

2.21 

3210 


1/160 

2.81 le-3 

1.79 

4885 

8.74 le-5 

2.06 

9390 

1.73 le-5 

2.07 

12620 

a 


a — 0.2 (Ay 

= As = h) 

a — 0.5 (Ay 

— A As — 

~ 10’ 5/ 

a = 0.8 (Ay 

_ h 

~ 50 

> As= A) 

Equal- 

1/10 

5.70 le-1 

- 

63 

2.86 le-2 

- 

65 

8.23 le-3 

- 

65 

height 

1/20 

2.01 le-1 

1.50 

208 

6.29 le-3 

2.19 

230 

1.62 le-3 

2.35 

230 


1/40 

4.04 le-2 

2.32 

700 

1.44 le-3 

2.12 

828 

3.49 le-4 

2.21 

860 

of (2.3) 

1/80 

9.59 le-3 

2.08 

2370 

3.43 le-4 

2.07 

2899 

8.03 le-5 

2.12 

3278 


1/160 

2.79 le-3 

1.78 

7887 

8.71 le-5 

1.98 

9752 

1.67 le-5 

2.26 

11521 

Equal- 

1/10 

5.81 le-1 

- 

35 

2.44 le-2 

- 

65 

7.77 le-3 

- 

69 

area 

1/20 

2.00 le-1 

1.54 

86 

6.38 le-3 

1.93 

182 

1.44 le-3 

2.43 

223 


1/40 

6.50 le-2 

1.62 

235 

2.82 le-3 

1.18 

529 

6.00 le-4 

1.26 

657 

of (2.3) 

1/80 

8.35 le-3 

2.96 

658 

3.84 le-4 

2.87 

1523 

1.52 le-4 

1.98 

1757 


1/160 

4.72 le-4 

4.15 

1903 

2.23 le-5 

4.11 

4112 

1.54 le-5 

3.31 

4315 

a 


a — 1.2 (Ai; 

' = As 

= lOh) 

a = 1.5 (Ay = As ■ 

= lOh) 

a = 1.8 (Ay = As = lO/i) 

Method 

h 

^max 

CO 

M 

^max 

CO 

M 

^max 

CO 

M 

Equal- 

1/10 

1.04 le-4 

- 

29 

9.00 le-5 

- 

29 

7.40 le-5 

- 

29 

height 

1/20 

4.46 le-6 

4.54 

59 

3.51 le-6 

4.68 

59 

2.61 le-6 

4.83 

59 


1/40 

1.88 le-7 

4.57 

119 

1.34 le-7 

4.71 

119 

8.95 le-8 

4.86 

119 


1/80 

7.86 le-9 

4.58 

239 

5.03 le-9 

4.73 

239 

3.04 le-9 

4.88 

239 


1/160 

3.26 le-10 

4.59 

479 

1.88 le-10 

4.74 

479 

1.02 le-10 

4.89 

479 

Equal- 

1/10 

1.04 le-4 

- 

29 

9.00 le-5 

- 

29 

7.40 le-5 

- 

31 

area 

1/20 

4.46 le-6 

4.54 

59 

3.51 le-6 

4.68 

59 

2.61 le-6 

4.83 

61 


1/40 

1.88 le-7 

4.57 

119 

1.34 le-7 

4.71 

119 

8.95 le-8 

4.86 

122 


1/80 

7.86 le-9 

4.58 

239 

5.03 le-9 

4.73 

240 

3.04 le-9 

4.88 

246 


1/160 

3.26 le-10 

4.59 

479 

1.88 le-10 

4.74 

482 

1.02 le-10 

4.89 

493 
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Table 2. The maximum errors {cmax), convergence rates (CO), and the total times of quadrature nodes 
being used (M) of Example at T = 5 and A = 1, by using the equal-height distribution as well as the 
equal-area distribution ones, for different a, h, Ay, and As, respectively. 


a 


a = 0.2 {Ay - 

--As = \) 

a = 0.5 {Ay - 

.As=|) 

a — 0.8 {Ay - 

= As=|) 

Method 

h 

^max 

CO 

M 

^max 

CO 

M 

^max 

CO 

M 

Equal- 

1/10 

4.36 le-2 

- 

458 

2.25 le-2 

- 

506 

5.60 le-2 

- 

571 

height 

1/20 

2.11 le -2 

1.05 

2137 

1.47 le-2 

0.61 

2472 

1.44 le-2 

1.96 

2790 


1/40 

9.64 le-3 

1.13 

10628 

7.16 le-3 

1.04 

11758 

7.85 le-3 

0.88 

12004 

of ( 2 . 2 ) 

1/80 

4.28 le-3 

1.17 

47104 

1.36 le-3 

2.40 

47721 

4.12 le-3 

0.93 

48101 


1/160 

1.87 le-3 

1.20 

189260 

5.77 le-4 

1.23 

191072 

2.12 le-3 

0.96 

192220 

Equal- 

1/10 

4.36 le-2 

- 

397 

2.40 le-2 

- 

454 

2.48 le-2 

- 

543 

area 

1/20 

2.11 le -2 

1.05 

1856 

1.55 le-2 

0.63 

2209 

1.44 le-2 

0.78 

2651 


1/40 

9.64 le-3 

1.13 

9070 

8.96 le-3 

0.79 

10651 

7.85 le-3 

0.88 

11513 

of ( 2 . 2 ) 

1/80 

4.28 le-3 

1.17 

41542 

1.67 le-3 

2.42 

43778 

4.12 le-3 

0.93 

46100 


1/160 

1.87 le-3 

1.20 

167412 

5.77 le-4 

1.54 

175082 

2.12 le-3 

0.96 

184222 

a 


a = 0.2 (Ay= 4, As = |) 

a = 0.5 {Ay= As = ^) 

a = 0.8 (Ay= As= A) 

Equal- 

1/10 

4.36 le-2 

- 

1159 

5.18 le-3 

- 

1243 

2.48 le-2 

- 

1181 

height 

1/20 

2.11 le -2 

1.05 

4193 

2.01 le-3 

1.36 

4538 

1.44 le-2 

0.78 

4173 


1/40 

9.64 le-3 

1.13 

19469 

1.18 le-3 

0.77 

16412 

7.85 le-3 

0.88 

14233 

of (2.3) 

1/80 

4.28 le-3 

1.17 

52550 

8.61 le-4 

0.45 

57446 

4.12 le-3 

0.93 

46425 


1/160 

1.87 le-3 

1.20 

180836 

5.77 le-4 

0.58 

196538 

2.12 le-3 

0.96 

142949 

Equal- 

1/10 

4.36 le-2 

- 

223 

5.14 le-3 

- 

471 

2.48 le-2 

- 

884 

area 

1/20 

2.11 le -2 

1.05 

585 

2.01 le-3 

1.36 

1471 

1.44 le-2 

0.78 

2921 


1/40 

9.64 le-3 

1.13 

1602 

1.42 le-3 

0.51 

4344 

7.85 le-3 

0.88 

8716 

of (2.3) 

1/80 

4.28 le-3 

1.17 

4394 

8.61 le-4 

0.71 

12096 

4.12 le-3 

0.93 

23458 


1/160 

1.87 le-3 

1.20 

12036 

5.77 le-4 

0.58 

32807 

2.12 le-3 

0.96 

57378 

a 


a = 1.2 {Ay = 

As = lOh) 

a = 1.5 {Ay = 

As = lOh) 

a — 1.8 {Ay = 

As = lOh) 

Method 

h 

67710 . 2 ; 

CO 

M 

^max 

CO 

M 

^max 

CO 

M 

Equal- 

1/10 

7.97 le-4 

- 

149 

2.82 le-3 

- 

149 

4.82 le-3 

- 

149 

height 

1/20 

2.44 le-4 

1.71 

299 

7.55 le-4 

1.90 

299 

1.27 le-3 

1.92 

299 


1/40 

6.66 le-5 

1.88 

599 

1.95 le-4 

1.95 

599 

3.27 le-4 

1.96 

599 


1/80 

1.73 le-5 

1.95 

1199 

4.95 le-5 

1.98 

1199 

8.27 le-5 

1.98 

1199 


1/160 

4.39 le -6 

1.98 

2399 

1.25 le-5 

1.99 

2399 

2.08 le-5 

1.99 

2399 

Equal- 

1/10 

7.97 le-4 

- 

156 

2.82 le-3 

- 

150 

4.82 le-3 

- 

155 

area 

1/20 

2.44 le-4 

1.71 

314 

7.56 le-4 

1.90 

303 

1.27 le-3 

1.92 

312 


1/40 

6.66 le-5 

1.88 

628 

1.95 le-4 

1.95 

609 

3.27 le-4 

1.96 

627 


1/80 

1.73 le-5 

1.95 

1256 

4.95 le-5 

1.98 

1221 

8.27 le-5 

1.98 

1256 


1/160 

4.39 le -6 

1.98 

2514 

1.25 le-5 

1.99 

2445 

2.08 le-5 

1.99 

2513 


19 



















Figure 3.1. The CUP time of the equal-height distribution as well as the equal-area distribution ones, for 
a = 0.5, T = 50, A = 1, h = 1/20, and Ay = As = h/2. 



Figure 3.2. The CUP time of the equal-height distribution as well as the equal-area distribution ones, for 
a = 1.5, T = 50, A = 1, h = 1/20, and Ay = As = lOh. 
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Table 3. The maximum errors (cmax), convergence rates (CO), and the total times of quadrature nodes 
being used (M) of Example at T = 4 and A = 1, by using the equal-height distribution as well as the 
equal-area distribution ones, for different a, h, Ay, and As, respectively. 


a 


a = 0.2 {Ay = As 

= lO/i) 

a = 0.5 {Ay = 

As = h) 

a = 0.8 (Ay = 

As = h) 

Method 

h 

^max 

CO 

M 

^max 

CO 

M 

^max 

CO 

M 

Equal- 

1/10 

1.08 le-1 

- 

122 

6.54 le-3 

- 

314 

1.07 le-2 

- 

331 

height 

1/20 

7.53 le-2 

0.52 

337 

3.58 le-3 

0.87 

1570 

5.24 le-3 

1.23 

1660 


1/40 

5.07 le-2 

0.57 

1404 

1.72 le-3 

1.06 

7189 

1.95 le-3 

1.43 

6877 

of (2.2) 

1/80 

3.31 le-2 

0.61 

7214 

9.00 le-4 

0.93 

28875 

5.49 le-4 

1.83 

27372 


1/160 

2.11 le-2 

0.65 

37504 

5.88 le-4 

0.61 

115585 

1.53 le-4 

1.84 

109480 

Equal- 

1/10 

1.08 le-1 

- 

119 

7.64 le-3 

- 

269 

1.07 le-2 

- 

295 

area 

1/20 

7.53 le-2 

0.52 

296 

4.05 le-3 

0.91 

1295 

5.46 le-3 

0.97 

1500 


1/40 

5.07 le-2 

0.57 

1073 

2.00 le-3 

1.02 

6162 

2.17 le-3 

1.33 

6221 

of (2.2) 

1/80 

3.31 le-2 

0.61 

5015 

9.00 le-4 

1.15 

24833 

6.12 le-4 

1.83 

24944 


1/160 

2.11 le-2 

0.65 

23817 

5.88 le-4 

0.61 

99550 

1.68 le-4 

1.86 

99922 

a 


q = 0.2 (Aj/= a. 

II 

<1 

Q = 0.5 (Ay = A 

, As= A) 

a = 0.8 (Ay = 

, As= A) 

Equal- 

1/10 

1.08 le-1 

- 

692 

4.05 le-3 

- 

314 

6.09 le-4 

- 

860 

height 

1/20 

7.53 le-2 

0.52 

2437 

1.80 le-3 

1.17 

1570 

2.81 le-4 

1.12 

3295 


1/40 

5.07 le-2 

0.57 

8468 

1.20 le-3 

0.59 

7189 

1.15 le-4 

1.29 

12458 

of (2.3) 

1/80 

3.31 le-2 

0.61 

2896 

9.00 le-4 

0.41 

28875 

4.37 le-5 

1.40 

45625 


1/160 

2.11 le-2 

0.65 

96830 

5.88 le-4 

0.61 

115585 

2.61 le-5 

0.74 

15926 

Equal- 

1/10 

1.08 le-1 

- 

125 

4.05 le-3 

- 

471 

6.09 le-4 

- 

707 

area 

1/20 

7.53 le-2 

0.52 

266 

1.80 le-3 

1.17 

1471 

2.81 le-4 

1.12 

2333 


1/40 

5.07 le-2 

0.57 

595 

1.20 le-3 

0.59 

4344 

1.15 le-4 

1.29 

6994 

of (2.3) 

1/80 

3.31 le-2 

0.61 

1378 

9.00 le-4 

0.41 

12096 

5.05 le-5 

1.19 

18771 


1/160 

2.11 le-2 

0.65 

3343 

5.88 le-4 

0.61 

32807 

2.38 le-5 

1.09 

45730 


Example 3. In this example, we examine the following initial value problem 

oDf’^x{t) = -x{t). (3.6) 

The initial values are given as a;(0) = 1 (and [e^*x{t)]'\^^^ = 0 if 1 < a < 2). It is known that J7^ the exact 
solution of this initial value problem is 

x{t) = e“^*£'a,i(-t“). 

Here the generalized Mittag-Leffler function is given by m 

oo ^ 

It is obvious that neither x(t) nor ((Df’^x^t) has a bounded first (second) derivative at t = 0 when 
0<a<l(l<a<2). Table shows the numerical accuracy of the equidistributing methods by choosing 
suitable Ay and As for different a G (0,1) when T = A and A = 1. And the equal-area distribution methods 
still show their advantages in computation cost. 


4 Conclusions 

For more deeply discussing the anomalous diffusion, sometimes the Levy waiting time distribution has to be 
tempered because of the finite lifetime of biological particles. Then the time derivative of the Fokker-Planck 
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equation describing this type of trapped dynamics is the tempered fractional operator. Since the tempered 
fractional operator is still nonlocal, generally the computation cost of numerically solving the tempered 
time-dependent problem is quadratically increasing with time t. This paper provides the predictor-corrector 
approach with theoretically proved convergence order min{l -|- 2a, 2} for the tempered fractional differential 
equation. The proposed predictor-corrector schemes on equidistributed meshes are detailedly discussed. In 
fact, the main contribution of this paper comes from the proposed equidistributing schemes which have 
linearly increasing computation cost with time t while keeping the accuracy; when a € (1, 2), the numerical 
results show that the convergence order can even exceed 4. And the larger the time is, the more benefits are 
obtained for the equidistributing methods. 


Acknowledgements 

The authors thank W.H. Deng for the discussions. This work was supported by the National Natural Science 
Foundation of China under Grant 11271173 and 11471150, the Fundamental Research Funds for the Central 
Universities (Northwest University for Nationalities) under Grant No. 31920130007, the Humanities and 
Social Sciences Youth Projects of the Ministry of Education of P R China under Grant No. 12YJCZH027 and 
No. 13YJCZH029, and the Humanities and Social Sciences Planning Projects of the Ministry of Education 
of P R China under Grant No. 11YJAZH053. 


References 

[1] B. Baeumera and M.M. Meerschaert, Tempered stable Levy motion and transient super-diffusion, J. 
Gomput. Appl. Math. 233 (2010) 2438-2448. 

[2] A. Cartea and D. del-Gastillo-Negrete, Fractional diffusion models of option prices in markets with 
jumps, Phys. A, 374 (2007) 749-763. 

[3] A. Gartea and D. del-Gastillo-Negrete, Fluid limit of the continuous-time random walk with general 
Levy jump distribution functions, Phys. Rev. E, 76 (2007) 041105. 

[4] W.H. Deng, Short memory principle and a predictor-corrector approach for fractional differential equa¬ 
tions, J. Gomput. Appl. Math., 206 (2007) 174-188. 

[5] W.H. Deng, Numerical algorithm for the time fractional Fokker-Planck equation, J. Gomput. Phys., 
227 (2007) 1510-1522. 

[6] W.H. Deng and C. Li, Numerical schemes for fractional ordinary differential equations. In: P. Miidla, 
(ed.) Numerical Modelling, pp. 355-374, InTech, Rijeka (2012), Chap. 16. 

[7] K. Diethelm and N.J. Ford, Analysis of fractional differential equations, J. Math. Anal. Appl., 265 
(2002) 229-248. 

[8] K. Diethelm, N.J. Ford, and A.D. Freed, A predictor-corrector approach for the numerical solution of 
fractional differential equations. Nonlinear Dynam., 29 (2002) 3-22. 

[9] K. Diethelm, N.J. Ford, and A.D. Freed, Detailed error analysis for a fractional Adams method, Numer. 
Algorithms, 36 (2004) 31-52. 

[10] N.J. Ford and A.C. Simpson, The numerical solution of fractional differential equations: Speed versus 
accuracy, Numer. Algorithms, 26 (2001) 333-346. 

[11] J. Gajda and M. Magdziarz, Fractional Fokker-Planck equation with tempered a-stable waiting times: 
Langevin picture and computer simulation, Phys. Rev. E, 82 (2010) 011117. 


22 



[12] A. Hanyga, Wave propagation in media with singular memory, Math. Comput. Model., 34 (2001) 1399- 
1421. 

[13] C. Li, W.H. Deng, and L.J. Zhao, Well-posedness and numerical algorithm for the tempered fractional 
ordinary differential equations, arXiv preprint arXiv:1501.00376, 2015. 

[14] M. Magdziarz and A. Weron, Competition between subdiffusion and Levy flights: A Monte Carlo 
approach, Phys. Rev. E, 75 (2007) 056702. 

[15] M.M. Meerschaert, Y. Zhang, and B. Baeumer, Tempered anomalous diffusion in heterogeneous systems, 
Geophys. Res. Lett., 35 (2008) L17403. 

[16] M.M. Meerschaert, F. Sabzikar, M.S. Phanikumar, and A. Zeleke, Tempered fractional time series model 
for turbulence in geophysical flows, Journal of Statistical Mechanics: Theory and Experiment 14 (2014) 
1742-5468. 

[17] 1. Podlubny, Fractional Differential Equations, Academic Press, New York (1999). 

[18] F. Sabzikar, M.M. Meerschaert, and J.H. Chen, Tempered fractional calculus, J. Comput. Phys., (2015), 
doi: 10.1016/j.jcp.2014.04.024. 

[19] S. Samko, A. Kilbas, and O. Marichev, Fractional Integrals and Derivatives: Theory and Applications, 
Gordon & Breach, London (1993). 

[20] M.G.W. Schmidt, F. Sagues, and I.M. Sokolov, Mesoscopic description of reactions for anomalous 
diffusion: a case study, J. Phys.: Gondens. Matter, 19 (2007) 065118. 

[21] A.B. White, On selection of equidistributing meshes for two-point boundary-value problems, SIAM. J. 
Numer. Anal., 16 (1979) 472-502. 

[22] L.J. Zhao and W.H. Deng, Jacobian-predictor-corrector approach for fractional differential equations, 
Adv. Comput. Math., 40 (2014) 137-165. 


23 


